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Abstract In this article, we present various numerical methods to solve multi-contact 
problems within the Non-Smooth Discrete Element Method. The techniques considered to 
solve the frictional unilateral conditions are based both on the bi-potential theory intro- 
duced by de Saxce et al. |2j and the Augmented Lagrangian theory introduced by Alart 
et al. [1]. Following the ideas of Feng et al. [3], a new Newton method is developed to 
improve these classical algorithms and numerical experiments are presented to show that 
these methods are faster than the previous ones and provides results with a better quality. 
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1 Introduction 

This is a first draft of a paper that will be submitted in a near future. 

The paper is organized as follow: in the next part, we present the equations to be solved 
for the Discrete Element Method, and the frictional contact law considered. In the third 
part, we first present two classical methods to numerically solve the full problem, the first 
one based on the bi-potential theory, and the second one on the Augmented Lagrangian 
theory. Then, we show how these methods can be enhanced using an appropriate Newton 
method. The last part on this article is devoted to the numerical experiments in order to 
show the main properties of these algorithms. 
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2 Problem Setting 



2.1 The equations of motion of a multi-contact system 

Classically (see for example |8l d [H]), the motion of a multi-contact system is described 
using a global generalized coordinate q (for Np particles, q G M'^^^p, where d = 6 for a 3D 
problem and d = 3 for a 2D problem). Due to the possible shocks between particles, the 
equations of motion has to be formulated in term of differential measure equation: 

Mdq + F*"*(f, q, q)dt = F^^*(t, q, q)dt + dR (1) 

where 

• M represents the generalized mass matrix; 

• F*"* and F*^^* represent the internal and external forces respectively; 

• dR is a non-negative real measure, representing the reaction forces and impulses 
between particles in contact. 

For the sake of simplicity and without lost of generality, only the external forces are 
considered in the following. The internal forces are neglected because the general case can 
be easily derived through a linearizing procedure. 

Then, for the numerics, the equation ([T| is integrated on each time interval [tk,tk+i], 
and approximated using a ^-method with 9 G]^, 1] for stability reason (see pHj). 
Therefore, the classical approximation of equation ([T| yields 

M(qfe+i - q„) = At{dFk+i + (1 - 9)Fk) + Rfc+i 

We will denote ql"'"'' = q^ +M"^ At(6'Ffc+i + (1 - 9)Fk) the free velocity (velocity when 
the contact forces vanish). Then, the first equation in ^ becomes 

qfc+i =qf"^ + M-iRfc+i. (3) 

In order to write the contact law, for a contact c between two particles (1 < c < Nc, 
where Nc is is the total number of contact), we define the local-global mapping 

u^ = P*(q,c)q 

R = P(q,c)r^ ^ ' 

where u*^ is the local relative velocity between the two bodies in contact and is the 
local contact forces (u'^, r*^ G M.'^ where d is the dimension of the problem, and P* is the 
transpose of matrix P). We also denote IP(q) the total-global mapping, for u and r in 
^dxNc (vectors composed of all relative velocity and contact forces respectively): 

u = P*(q)q 

R = P(q)r ' 
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In the discretization, a prediction of q is computed to estimate the mapping P(q) (see 



equations (18) and (19) in the following). 



Using the equations ^ and the discretization of the motion of a multi-contact 
system, with frictional contact between particles can be written: 

lawc{ul^i,rl^^) = .true. Vc G {1, 2, A^J 
where W = P*M^-'^P is the Delassus operator, and u^*"^*^ = P*q^''^^ is the relative free 



velocity. Notice that a Newton impact law is also considered (see [TT\ and equation (20) 
in the following), that modify u^, and ujj''^'^ by and u{^^^ respectively. 

The second equation in ([g]) is the implicit frictional contact law that is in our case the 
classical Signorini condition and Coulomb's friction law. 

2.2 The frictional contact law 

In the local coordinates system defined by the local normal vector n and the tangential 
vector t _L n, any element u and r can be uniquely decomposed as u = u^n + and 
r = r„n + rj respectively. In these coordinates, the unilateral contact law can be stated 
using the Signorini's conditions (see figure [T] for a graphical representation): 

Un > 0, r„ > 0, ii„r„ = 0. (7) 



Contact 



No contact 



Figure 1: The Signorini conditions 



On the other hand, the Coulomb's law of friction can be stated using the algorithmic 
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form (see figure [2] for a graphical representation): 
If r„ = then u„ > 

Else if > and \\rt\\ < fiVn then u = 



! No contact 
! Sticking 



Else r„ > and ||r(|| = /ur„ then 3A > such that = ! Sliding 



(8) 



Sliding 



Sticking 



Figure 2: The Coulomb conditions 



For a given friction coefficient fi, let Kf^ be the isotropic Coulomb's cone, which defines 
the set of admissible forces (see figure |3]): 



= {r = r„n + rt : \\rt\\ - /ir„ < 0} 



(9) 













^ 



Figure 3: The Coulomb's cone 
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The previous law can be also written: 



r If r„ = 



then Un > 



! No contact 



Else if r G I{K^) 



then u = 



! Sticking 



(10) 



Else r„ > and r G B{K^) then 3A > such that = ! Shding 

where I{K^) and B{K^) are respectively the interior and the boundary of the cone K^. 

3 Numerical Resolution of the contact/friction problems 

We will describe in this section the numerical algorithms that will be considered in the 
following. Generally, to solve the problem ([6|, the numerical algorithms considered are 
based on two levels: the global level where the equations of motion are solved, and the 
local level devoted to the resolution of the contact law. 

3.1 Resolution of the global problem : the Non Linear Gauss Seidel 



In this paragraph, we describe the algorithm used at the global level to solve the problem 
([6|. Following the ideas of Jean and Moreau [HI [iTj, we use the non-linear Gauss-Seidel 
algorithm which is the most commonly used. It consists in considering successively each 
contact until the convergence. The numerical criterion used to state the convergence will 
be studied latter in the paper. 

This method is intrinsically sequential but it is possible to used a simple multi-threading 
technique which consists in splitting the contact loop into several threads. This method 
has been studied in |I6| in the case where the local algorithm is based on the Augmented 
Lagrangian method. 

Notice that it is also possible to consider at this stage more sophisticated method such 
as a conjugate gradient type method (see for example |14j). 

3.2 The standard bi-potential based method (SBP) 

In this paragraph, we provide a first method to solve the contact problem, at the local level 
(contact point between two particles). The method is based on the notion of bi-potential, 
introduced by de Saxce et al. |2]- 

Using the bi-potential framework, it can be shown (see for example O IH [6l ET]) that 
a couple (u, r) verifies the Signorini-Coulomb contact rules if 



Method (NLGS) 



6c(v, s) -I- V • s > 6c(u, r) -I- u • r = 



Vv,s 



(11) 
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where be is the bi-potential 

6c(-u,r) = ^K+{un) + ^K^ir) + fJ-rn\\ut\\ (12) 

and stands for the indicatrix function of the set C: ^c{x) = if x € C, ^c{x) = 
+00 if X ^ C. 

Consequently, the contact law can be written in a compact form of an implicit subnor- 
mality rule (or a differential inclusion rule): 

-u£dM-u,r). (13) 

Then, for a contact c, at a NLGS iteration i, knowing the relative velocity u^'*, the 
algorithm to compute r^'*"*"^ from r^'* is based on the minimization of the bi-potential (see 
for exemple [3], page 51), using the inequality: 

6c(-u''\ r) + u""'' • r > bc{-u'''\r^''+^) + u'''' • r^'^+^ Vr G (14) 
or g(r) > g{r'^'^~^^), Vr S K^, if we denote 

9{r) = *M+«'*) + ^A' (r) + f^rJurW + Ci^'* • r. (15) 



The minimization of (14) is classically realized using a projected gradient projection 
(Uzawa method) without considering the singular term ^^+(un^). This minimization can 
also be viewed as the proximal point of the augmented force r — pu, with respect to the 
function r 1— t- pbc{—u, r) (see for example [21 HI |6]): 

r = prox{r — pu, pbd—u, r)). 



More precisely, the Uzawa method leads to compute the augmented force r'^'*^^ = 
r'^'* — pV5(r^'*), where g is the differential part of g: 

y-gir"'') = Vr(Mr„||u^''|| + ii'^'^ • r) = ;u||a^'*||n + 

and to consider the force at next step as a projection of the augmented force onto the set 
of admissible force r'^'*"'"-'^ = proj{T'^'^~^^ , K^j^), that provides equations (21) and (22) in the 
resolution algorithm of the global problem. The proj{T^''^~^^ , K^) stands for the orthogonal 
projection over the convex K^^ that can be computed exactly (see j^). 

This algorithm will be referred as the SBP (Standard Bi-Potential) method above and 
throughout. 

For a sake of simplicity, we denote hereafter the descent direction 

D'^'* = //||Ui'*||n + ii^'\ 
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Remark 1 A first improvement of this method could be to compute the optimal step p^'^ . 
To do so, we have to minimize 

p^5(r'^'^-pD'='^), (16) 

or, more precisely, 

p ^ ^M+{un') + ^K.ir''' - pD^'O + pir'rl' - pB''' ■ n)\\u'/\\ + a^-' • (r^-' - pD^'O 
= ^M+(nS') + - pD^'O - pB''' • {p\\u/\\n + u^'') + Cte 

= ^^+{u'k') + ^K,{r^^' - pD'^'O - pIlD^'^f + Cte. 

(17) 

We can observe that this method do not permit to choose an optimal parameter p since 
g, as a function of p, is linear, excepted in the case where D"^'* ^ K^. A solution could 
be to modify the function g, for example by replacing u^'* by a prediction of u'^'*+i using 
the equations of the dynamics. Unfortunately, this method do not provides good numerical 
results. 



Then, the standard bi-potential based algorithm (SBP) can be written (see |21] for 
example) : 

• Loop on the step time k 

— Prediction of a position (for the computation of the local-global mapping): 

At. 

q^+i = qfc + ^qfc; (18) 

— Initialization of the motion: q^l_^_i = q{'~'^'^ (initialization of the contact forces 
with R = 0). 

— Loop on i > (NLGS), until convergence 

* Loop on the contacts c: 

• Computation of the local-global mapping 

= P*(q,+ i,c)qfc ; u'='+^ = P*(q,+ i , c)qi,+i (19) 

• Newton shock law 

c,+i . — c,+i . — 

~c,i _ + e-nUn . ~c,i _ + 671 Ut /^r,\ 



Prediction of the reaction: 



+ {U'^n' + l^WK'W)^ (21) 
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• Correction of the reaction: 

j.c,i+i ^ proj{T'''^+\K^) (22) 

• Actualization of the generahzed displacement: 

^kX\ = + M-^E ^(q^+i ' + E ^(q^+i ' (23) 

a<c a>c 

* End of the loop on contacts c. 

— End of the loop on i of NLGS when the convergence is reached: q^+i = 

— Actualization of the generalized displacements: q^+i = q^+i + ^qfc+i 

• End of the loop on the step time k. 

Remark 2 Notice that only one iteration of the Uzawa algorithm at the local level is 
considered. Various previous studies (see for example 191) show that there is no significant 
improvement of the method if several iterations of the Uzawa algorithm are considered at 
this stage. 

3.3 Newton method and enhanced bi-potential method (EBP) 

We introduce in this section a Newton method in order to speed up the convergence of the 
computation of the solution. This method has been already used, especially in the case of 
the augmented lagrangian method developed by Alart et al. [l], and the ideas presented 
in this article follows those of Feng et al. p] and have been adapted to the problem of 
the discrete element method. The main idea of this technique is to find the solution of 
the optimization problem, not as a minimum of a functional, but rather as a zero of a 
function, using the Euler equation of the problem. Then a standard Newton method can 
be developed to solve this Euler equation. 

The technique is first described in the case of the bi-potential framework, and will 
adapted to the augmented lagrangian method farther. 

We recall that the local problem that has to be solved, for each contact c can be written 

^k+i-^k +2^ Wear Vc=l, A^e (24) 

, = proj{T'',K^) 



where = — p(/x||u^||n + u) is the augmented reaction (see 21), and 
P*{^.k+'^ c)^^^-^(^A:+i a) local Delassus operator. 
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This problem can be written equivalently 



U 



fc+1 



~ c,free 



U 



^ t^ear" = 



a=l 



Vc= 1, iVe 



(25) 



^ - proj{T'',K^) = 



Reminding now that we want to use a Newton algorithm to solve theses equations 
inside the Non Linear Gauss Seidel loop on the variable i, we define now, for each contact 
c = 1, ...jNc, the function 



fcix) 



I 



V 



U^'* - u 



c,free 



a=l 



\ 



J 



where : 

• the vector Z'^ is the error on the prediction of the reaction 

Z^'*(r^'\ii 



(26) 



Xc = (r-^'SCi^'*)*, 



X = {xi,X2, ■■■,XN,Y 



Remark 3 The first equality in the relation fix) = is the equation of motion for the 
bodies in contact, and the second relation is the frictional Coulomb law between the bodies 
in contact, written within the bipotential framework. 

Then we have to write a Newton algorithm to solve the problem f{x) = 0. This 



algorithm can be written, for a contact c, by substituting equations (21) and (22) in 
algorithm (SBP) by the followings: 



• Initialization: 



X 



v° = u'='') , 1 = 



Loop on £, until convergence: 

— = — /)(/z||v|||n + 

— Resolution: 



dx' 



ix' 



Axc = -fc{x' 



(27) 
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— Actualization: Xc^^ = Xc + ^Xc 
End of the loop on £ until convergence, u'^''"*"-'^ = and r'^''^~^^ 



r . 



Remark 4 This algorithm needs more than one iteration at each Non Linear Gauss Seidel 
iteration to he efficient. As a consequence and compared to the Uzawa algorithm, the 
solution in the Newton algorithm is controlled by both the local (iteration £) and global 
convergence criteria (iteration i, see 



The local convergence criterion for the Newton algorithm is defined by: 
e%eUXi) = h'- '■^^ - Wr% + \\/-proj{r\K,)\\ 
This criterion measure fdxi) ^^^^ has to be sufficiently small. 



(28) 



The matrix 



represents the tangential matrix of the local equations for the 
contact c. This matrix is of dimension 6x6 for a 3 dimensional problem, and 3x3 for a 2 
dimensional problem. For a 3 dimensional problem, the general form of this matrix is the 
following: 



where 



drn 



dfc 

dXc 



dZ, 



ix) 



dZ, 



-W Id 



3x3 



drt 



dr 



t2 



Br 



Br 



dZ, 

dVn 



dZr 



dvt^ 



dZr 

dv 



t2 



(29) 
(30) 



The matrices Ac and Be takes different forms according to the contact status: 



First case: sliding contact. 
In that case, we have 



then 



and 



Ai||Tt|| > -Tn 

Proj{T,K^) = T - 
= = /o(/"l|vt^||n + v^) + 



In II > fJ'Tn 



1 + 



n 



/in 



n 



fin 



The computation of the derivatives of Zc provides the matrices Ac and Be'- 



dZc _ fi / _Tt_ 

drn 1 + Vikil, 
dZr _ Tt, / jn_ 
drt, ~ {l + fi^)\\Tt\\ Vllnl 



/in 



fin + 



mil - 
l+/i2 



^1 
Inll' 
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^ f Tt \ llnll - l^rn / t2 Tta ^ 

drt, il + iJ,^)\\Tt\\\\\Tt\\ ^"y'^ 1 + Vllrtll Urtf""*, 
dZc , pfi f n 

— pn+ -— y J r. - jJLU 



dvn 1 + win 



For a 2D problem, these computations yields: 
dZc n 



dvn 1 + /x^ 
dZc 1 



(/xn — dr't) 
{—nOrn + 1) 



5rt (1 + fi^) 
dZc pp, 



dvt 1 + p^ 
where 6^ = sign(vt) and Oj. = sign(rt). 



• Second case: sticking contact. 
In that case, we have 

Mllnll > -Tn llnll < pTn 

then 

Z, = p(M||vf||n + v'=) 
and the computation of the derivatives of reads: 





03x3 




dZ, 






dVn 


= pn 




dZc 






dvt^ 


= Pl^- 




dZ, 




Ut2 


dvt2 


= Pl^- 





n + pti 
n + pt2 
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For a 2D problem, these computations leads to: 



dVn 

dvi 



2x2 



pfiOi.n + pt 



• Third case: no contact. 

In that case, the matrices Ac = Idsxs and vanishes, and xP'^ = 

3.4 Resolution of the linear system 

Generally, the drawback of a Newton is the computational cost of the linear system to be 
solved at each iteration. Here, the particular form of the tangent matrix allows the use of 
a condensation technique. More precisely, the linear system to be solved can be written: 



-W 7^3X3 



5r 
(5v 



(31) 



The first equation yields 6v = — f + W6r, and introducing this equality is the second 
equation leads to solve the linear system 



{Ac + BcW)5r = -g + BJ. 
This properties halves the size of the linear system to be solved. 



(32) 



Remark 5 A drawback of the bi-potential framework is that, due to is specificity, it is 
rather difficult to consider fully coupled problems, where the contact law and another phe- 
nomena, such as electricity or thermic effects are strongly coupled. The other method 
presented in this paper has a better property from this point of view because it is based on 
a more standard mathematical background in the theory of optimization. 



3.5 Newton method and enhanced augmented lagrangian method, (SAL) 
and (EAL) 

In [1], Alart et al. propose another method to solve the frictional contact problem. This 
method has been also used with various improvement (parallelization, conjugate gradient 
method for example) to solve multi-contact problems |14| [TSl [TBI \T7\ [TS] . Even if the 
coupled frictional contact problem is not an optimization problem anymore, it is always 
possible to formally formulate a "quasi"- optimization problem, for which the constraint 
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set depends on the normal components of the solution as a parameter. The solution is 
then searched as a saddle point of a "quasi" augmented Lagrangian of the problem. 

More precisely, the global problem on all unknowns that has to be solved at each time 
step (in place of equation (24)) has the following form: 

free 



U = 



+ 



(33) 



r > 0, u > 0, r • u = 0. 



In order to solve this problem, for a given r G M^^^'^, one can define the cartesian 
product of infinite half cylinder with section equal to the ball B{0, ^r^) of radius /xrc by: 

Nc 



c=l 



and then, the granular type frictional contact problem is given by 



1 



free 



(34) 



r G argmmrgc(Air) - r • Wr + u-* 

and the projected gradient method the minimize this problem reads (for each iteration i 
of the NLGS algorithm): 

r*+i = projir' - p(u^^'^'= + Wr*), C{nr'+^)), (35) 

or r*+^ = proj{T^~^^, C(/ir*+-^)), with t^~^^ = r* — pu\ u* = u-'^^^^ + Wr*. This algorithm 
will be referred to hereinafter as the SAL (Simple Augmented Lagrangian) method. 

Notice that this method is very closed to the SBP method. More precisely, for a contact 
c, only the descent direction u'^'' + ;u||u^'*||n in (21) is replaced by u'^'* and the projection 
],c,i+i = p'roj(T'^''^^^ , Kfj_) in (23) is replaced by 



/ c,i+l 



c,i+l 



max{0, r, 



c,i+l\ 
n ) 



c,i + l 

'j. 

C,! + l| 



Remark 6 On the contrary, it is possible to see the algorithm developed from the bi- 
potential formalism as a slight modification of the algorithm above. Indeed, it is only 

necessary to change the set C{r) by K, = W^^, and to change the descent direction u'^'* 

c=l 

by u'^'' + /x||u^'' ||n which remains a descent direction for the SAL method, since 
VJ(r^'*+i) • D^'^ = -||u"i|2 - li"'* • (/x||uti|)n = - ^u^^'W^/W 

which is negative since pi G [0, 1]. 
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Then, acting by analogy, we can develop a Newton method to find the minimum of J 
by seeking the solution as a zero of the function /(x) where, for a contact c 



/c(x) 



V 



a=l 



\ 



J 



the vector 7f is the error on the prediction of the reaction 

Z%v'^, u^+i) = - proM^MfiT^^,)), 



(36) 



and the set Cc(/xr'^) is the set of admissible forces Cc(/xr"^) = M"^ X B{0, Vc). This method 
will be refered as the EAL (Enhanced Augmented Lagrangian) method hereafter. 



Then, as bellow, we have three cases in the computation of the tangent matrix 

• First case: sliding contact (r„ > 0, tj > /ir„) 

We have: proj{T'',Ccipr'')) = r„n + ||^|Ur„t and Z,. = pi;„n - ^^t„ + r^. 

The computation of the derivatives of Zc provides the matrices Ac and Be- 

n 



dZc 
drn 
dZc 

drti 

dZc 
drt2 
dZc 

dVn 

dvt^ 
dZ, 



t2 



~I^U „ 
ll'^tll 

ti - fJ,Tn 
t2 - fJ-Tn 
n + fi 

-PIJ'Tn 



_t2_ 
In I 



^1 
Inll' 

^2 
Inll- 



■n 



n 



n 

WnW 
h_ 

\\n\\ 

_t2_ 



^1 



For a two dimensional problem, these computations yields 





—flOr 1 



Br 



1 

flOr 



• Second case: sticking contact (r„ > 0, tj < /7,r„) 
proj {t^ , Ccil^T'^)) = and the computation of the derivatives of Zg reads 
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- Ac = 03x3 

- Bc = /oldaxs 



Third case: no contact (t„ < 0) 

proj{T^,Cc{fiT'^)) = 0, then the matrices Ac = Idsxs and Be vanishes, and Xc^^ 


3.6 The global stopping (convergence) criterion 

We present in this paragraph the convergence criterion on the global non linear Gauss- 
Seidel iterations. This criterion, developed from that proposed in has been extended 
in the case of the Newton and bi-potential (EBP) method, where some term are naturally 
vanishing in the original Uzawa and bi-potential (SBP) method. This criterion Egiob has 
been written in such a way that if the solution verify that Egiob is sufficiently small, then 
this solution has good properties on the equation of motion and Signorini Coulomb contact 
law. Consequently, this criterion stays valid for the methods developed with the augmented 
lagrangian (SAL and EAL methods). 

This criterion can be stated: 

1 

^glob — ~jrr ^ ^ [^motion ~^ ^proj ~^ ^bc ~^ ^pen\ i^'^) 



c=l 



where: 



^motion = W^" - "mil where = u^'* + Yfl=i Wear", so Emotion measures the error 
on the equation of motion (see equation ( |24[ ), this term vanishes for the SBP and 
SAL method); 

^proj ~ \/\\^^ ~ Pfoj{r'^, -^A')P is the error for the projection on the Coulomb cone 
(vanishing for the SBP method); 



u^-r^-|-/ur^||uj 



is the absolute value of the bi-potential that has to vanish if 



and only if the couple (u"^, r'^) verifies the Signorini Coulomb contact law (see formula 



111; 



^ven — ~ min(0, vf^) is the value of the penetration. 



Remark 7 One can notice that is absolutely necessary to verify in the criterion that there 
is no penetration, because nothing in the presented algorithm ensures that is condition is 
verify at the end of the loop. Moreover, if this condition is not satisfied, the rest of bi- 
potential can be negative or equal to zero, even if the couple (u, r) is not a solution. 
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4 Numerical results 



We present in the section three numerical examples with an increasing complexity. 

In these computations, the descent parameter p is taken in such a way that the result 
is optimal, in terms of time computing. Denoting p = for the SBP and the 

SAL methods, we have chosen p = 0.6p, whereas for the EBP and the EAL methods, it 
is better the take p = p. We recall that it has been show that, for the the bi-potential 
method (see for example [3]) and the augmented lagrangian method (see for example jl4J), 
the parameter p has to verify p < 2p in order to ensure the convergence. Generally, for 
these two methods, the convergence is very sensitive on this parameter. We will show in 
the last paragraph of this study that for the EBP method, the parameter p can be taken 
in a large range around the value p without changing dramatically the convergence of the 
method. 

At each iteration of the NLGS algorithm, the Newton algorithm is stopped either if 
the convergence is obtained (e^ewt — the number of iteration of the Newton 

algorithm reached 100 when there is no convergence. 

4.1 Ball sliding on a plane 

In this first example, we consider a ball placed on a table with an initial horizontal velocity 
equal to 1.5 m- s~^. The ray of the ball is equal to 5 • 10"'^ m, and the friction coefficient 
between wall and ball is equal to p = 0.7. The time step of discretization is equal to 10~^ 
s. In this experiment, the ball ffist slides on the table, and then the ball rolls without 
sliding. The global sloping criterion is equal to e^/ofe = 10"^*^. 



Figure 4: Example 1 - A ball is launched with an initial horizontal velocity (left). First, 
the ball slides. Then, the ball rolls without slipping (right). 
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Method 


Number of 
NLGS iterations 
(last time step) 


Error 

^glob 

(last time step) 


Total 
CPU time (s) 


SBP 


18 


0.384 • 10"^" 


9.44 


SAL 


18 


0.384 • 10-1° 


9.28 


EBP 


1 





8.83 


EAL 


1 


0.175 • 10-13 


8.78 



Table 1: Comparison of the results obtained by the four methods on the first example 
(after the 2000th time step). 



SBP 

0'' I ' ' ' ' ' 




10 4 ' ' ' ' — ' — ' — ' — ' — 4- 

10 10 
Non linear Gauss-Seidel iterations 



Figure 5: Example 1 - Convergence for the standard bi-potential based method, 5* iter- 
ation 
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SAL 

10'^ I ' ' ' ' ' 




1 o""^*!- ' ' ' ' ' ' ' ' L_ 

10 10 
Non linear Gauss-Seidel iterations 



Figure 6: Example 1 - Convergence for the standard augmented lagrangian method, 5 
iteration 




Figure 7: Example 1 - Convergence for the Newton and bi-potential method, 5 iteration 
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We can observe from these numerical results that the error coming from the projection 
is very small for the four methods. The Standard Bi-Potential (SBP) method and the 
Standard Augmented Lagrangian method (SAL) give very closed results, both in term of 
quality (see figures [s] and |6]) and in term of time computing (see table [T]). Nevertheless, we 
can notice that the time computing is smaller with the SAL method, because there is less 
computations at each iteration (no term such as ||u(|| and projection easier to compute 
for example). The Enhanced Bi-Potential method provides better results, both in term 
of quality (see figure [?]) and in term of time computing (6.5% better). The Enhanced 
Augmented Lagrangian method converges after the first Non Linear Gauss Seidel iteration 
for every time steps, and consequently, this is the faster method on this example (7% faster 
than the SBP method). 



4.2 Sedimentation of 4 balls in a box 



In this second experiment, we consider the sedimentation of 4 balls of radius ranging from 
4 • 10~^ m to 5 • 10~^ m. For the computations, the time step of discretization is equal 
to At = 10~^ s., and the Non linear Gauss-Seidel loop is stopped either if the the global 
stopping criterion on the NLGS method is equal to Sgiob = 10"^'^, or after 5000 iterations if 
there is no convergence (this case never occurs in this experiment). The friction coefficient 
between the balls and between the balls and the walls is equal to /i = 0.3. 




Figure 8: Example 2 - Sedimentation of four balls under the gravity effect. 



Like in the previous simulation,, the SBP and the SAL method methods provide very 
similar results (see figures [9] and 10). For these two methods, we can notice that here the 
global error is essentially due to the penetrations. The SAL method is 2% faster than the 
SBP method (table [2]). 



Results obtained by the EBP method are better (figure 11), and here the overall error 
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Method 


Number of 


Error 


Maximal 


Total 




NLGS iterations 


^glob 


penetration 


CPU time (s) 




(last time step) 


(last time step) 


(last time step) 




SBP 


305 


0.949 • 10"^^ 


0.310-10-11 


2.92 


SAL 


301 


0.980 • 10-12 


0.340 • 10-11 


2.87 


EBP 


161 


0.635 • 10-12 


0.641 • 10-12 


2.59 


EAL 


158 


0.973 • 10-12 


0.208 • 10-1^ 


2.43 



Table 2: Comparison of the results obtained by the four methods on the second example 
(after the 1000th time step) 



Bi-potential 

Penetration 




Non linear Gauss-Seidel iteration; 



Figure 9: Example 2 - Convergence of the non-linear Gauss-Seidel iterations for the stan- 
dard bi-potential based method (1000th time step) The two last curves overlaps, showing 
that the global error is governed by the error of penetration. 
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SAL 




1 L ' ' ' ' ' — ' ' ' I ' ' ' ' ' — ' ' ' I ' ' 1 

10 10 10 

Non linear Gauss-Seidel iterations 

Figure 10: Example 2 - Convergence of the non-linear Gauss-Seidel iterations for the 
standard augmented lagrangian method (1000th time step). The two last curves overlaps, 
showing that the global error is governed by the error of penetration. 



EBP 

Bi-potential 




Non linear Gauss-Seidel iterations 

Figure 11: Example 2 - Convergence of the non-linear Gauss-Seidel iterations for the 
Newton and bi-potential method (1000th time step). The two last curves overlaps, that 
shows that the global error is governed by the error on the equations of motion. 
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EAL 



— Bi-potential 




ition 

ration 

motion 



Non iinear Gauss-Seidei iterations 



Figure 12: Example 2 - Convergence of the non-linear Gauss-Seidel iterations for the 
Newton and Augmented Lagrangian method (1000th time step). The two last curves 
overlaps, and the other ones does not appear on the figure because the corresponding 
errors are lower than 10~^^. 



is governed by the error on the equations of the motion. The EBP method is 11% faster 
than the SBP method, and the penetration is 5 times smaller. In this example the EAL is 
the faster method (16,8% faster than the SBP method), and the penetration is very small 



(see figure 12). 



4.3 Sedimentation of 500 balls 



In this example, we consider the sedimentation of 500 balls (see figure 13) of radii ranging 
from 2.5 • 10~^ m to 5 • 10~^ m, the time step of discretization is equal to At = 5 • 10~^ s, and 
the Non linear Gauss-Seidel loop is stopped if the global estimator (37) verifies Egiob < 10~^^ 



or after after 5000 iterations if there is no convergence. The friction coefficient between 
the balls and between the balls and the walls is equal to fj, = 0.3. 
The results in table [3] are obtained after 1000 time steps. 

In this example, the difference between methods SBP and SAL on the one hand, and 
the method EBP and EAL on the other hand is larger (see table |3]). We can notice that the 
SAL method is 10.95% faster than the SBP method, and the EBP is 21.83% faster than 
the SBP method. Here, the EAL method is no longer the faster one, but the penetration 
is very small. Again, for the two first methods the global error is essentially due to the 
penetrations whereas the two last methods, the error is essentially due to the failure to 
follow precisely the equations of motion. 
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Figure 13: Example 3 - Zoom on balls falling under the gravity effect. Initial configuration 
on the left, final configuration on the right. 



Method 


Number of 


Error 


Maximal 


Total 




NLGS iterations 


^glob 


penetration 


CPU time (s) 




(last time step) 


(last time step) 


(last time step) 




SBP 


5000 


0.119-10-6 


0.213-10-5 


1092.95 


SAL 


5000 


0.135 • 10-^ 


0.533-10-5 


973.31 


EBP 


5000 


0.156-10-6 


0.286-10-6 


854.31 


EAL 


5000 


0.101 - 10-6 


0.390 - 10"^^ 


916.65 



Table 3: Comparison of the results obtained by the four methods on the third example 
(after the 1000th iteration, Nmax = 5000 iterations) 
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SBP 




10''^- 



lO''"- 



10° 10' 10^ 10^ 

Non linear Gauss-Seidel iterations 



Figure 14: Example 3 - Convergence of the non-linear Gauss-Seidel iterations for the 
standard bi-potential based method (1000th time step). The two last curves collapse. 



SAL 




0'"- 



10° 10' 10^ 10° 

Non linear Gauss-Seidel iterations 



Figure 15: Example 3 - Convergence of the non-linear Gauss-Seidel iterations for the 
standard augmented lagrangian method (1000th time step). The two last curves collapse. 
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Figure 16: Example 3 - Convergence of the non-linear Gauss-Seidel iterations for the 
Newton and bi-potential method (1000th time step). The two last curves collapse. 




Figure 17: Example 3 - Convergence of the non-linear Gauss-Seidel iterations for the New- 
ton and Augmented Lagragian method (1000th time step). The two last curves collapse. 
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4.4 Discussion on the descent parameter p 

We consider again the third example solved by the Newton and bi-potential method {etot = 
10~^, maximal number of iterations of Newton method equal to 100, £Newt = 10~^, 500*^* 
time step). Here, we take p = and we consider p = ap, for various values of a. 



a 


Number of NLGS 


Maximal 


Total CPU 




iterations 


penetration 


time (s) 




(last time step) 


(last time step) 




5 


652 


0.110 • 10"^ 


65.08 


2 


414 


0.177-10-6 


55.86 


1 


750 


0.149 • 10-6 


53.95 


1 

! 

5 


812 


0.634 • 10-6 


78.43 


667 


0.219 • 10-5 


176.14 



Table 4: Comparison of the results obtained for various values oi p = ap on the third 
example (after the 500th iteration, N^ax = 5000 iterations) 

These results show one of the main advantage of the EBP method. Indeed, one can 
notice that in table |4j the CPU time and the quality of the solution are very similar if a 
is equal to 1 or 2. Even if a is equal to five, the convergence is not to damaged. In that 
case, one remain the the SBP and the SAL are no longer convergent. If the parameter 
a is small, the method converges but the convergence rate is very small. One can notice 
that the EAL method is much more sensitive about the parameter a, essentially in the 
convergence of the Newton method. 

5 Conclusion 

The results presented show that, using an appropriate Newton method, it is possible to 
improve the computational time over that 20% compared to the standard methods. More- 
over, one principal drawback of that type of methods, that is the dependance of the results 
on the parameter p does not exist anymore. 

In the future, this method will be extended to the case of a contact law with adhesion. 
This improvement will be realized in a near future. 
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